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Abstract 

We consider stochastic descriptions of chemical reaction networks in which there are both fast and slow 
reactions, and for which the time scales are widely separated. We develop a computational algorithm that 
produces the generator of the full chemical master equation for arbitrary systems, and show how to obtain a 
reduced equation that governs the evolution on the slow time scale. This is done by applying a state space 
decomposition to the full equation that leads to the reduced dynamics in terms of certain projections and the 
invariant distributions of the fast system. The rates or propensities of the reduced system are shown to be the 
rates of the slow reactions conditioned on the expectations of fast steps. We also show that the generator of 
the reduced system is a Markov generator, and we present an efficient stochastic simulation algorithm for the 
slow time scale dynamics. We illustrate the numerical accuracy of the approximation by simulating several 
examples. Graph-theoretic techniques are used throughout to describe the structure of the reaction network 
and the state-space transitions accessible under the dynamics. 

Keywords'. Stochastic dynamics, reaction networks, graph theory, singular perturbation 


1 Introduction and background 

While singular perturbations techniques and the quasi-steady-state approximation (QSSA) have a long history 
of use in deterministic descriptions of chemical reaction kinetics (cf. [Lee & Othmer, 200^ for a review), it 
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was apparently not applied to discrete stochastic descriptions of chemical kinetics until Janssen proposed a 
method for adiabatic elimination of fast variables in stochastic chemical reaction networks [Janssen, 1989a[ 
Janssen, 1989b| |. Janssen began with a master equation description and, using projection techniques related to 


volume expansions, obtained a reduced master equation and showed in some examples that an intermediate 
chemical species can be eliminated in a network if a reaction occurs much faster than the others. Moreover, he 
also gave examples in which the master equation cannot be reduced via these techniques [Janssen, 1989b|| . It is 
sometimes assumed that the results of a deterministic reduction of a reaction network yields correct results for 
the stochastic description, but examples show that this is not correct [Thomas et al., 20111 . 

The Gillespie algorithm [Gillespie, 2007] is the most widely-used algorithm for simulating stochastic re¬ 
actions but it can be very inefficient when there are multiple time scales in the reaction dynamics. Stochastic 
simulation algorithms of multiscale reaction networks were developed more or less simultaneously by Haseltine 
and Rawlings ( [Haseltine & Rawlings, 2002| ) and Rao and Arkin ( [Rao & Arkin, 20031 ). The former authors 
formulated the changes in the species numbers n in terms of extents r] by defining n{t) = n(0) -h and 

divided the reaction extents into those of slow reactions and fast reactions. Rao and Arkin ( [Rao & Arkin, 20031 ) 
divided the set of species into ‘primary’ {y) and ‘intermediate or ephemeral’ (z), and assumed two conditions 
on the conditional variable z\y, (i) the conditional variable z\y is Markovian and (ii) z\y is at quasi-steady- 
state on the slow time scale. However, the first assumption was not justified and one cannot classify species 
as slow and fast - rather it is reactions that are slow or fast and species can participate in both. Mastny et al. 
( [Mastny et al, 2007| ) applied singular perturbation analysis to remove the QSS species for a number of model 
systems having small populations, and for networks where non-QSS species have large populations, they also 
utilized the fi-expansion to reduce the master equation. More recent work deals with the role played by the 
reduction on the level of noise in the solution [Srivastava et aL,20llt 

Several others have proposed and implemented variations of these two approaches for hybrid simulations of 
stochastic systems. These include partitioning the system dynamically [Sails & Kaznessis, 2005[ and solving 
the fast reactions using an SDE approximation. The effective reaction rate expressions for Hill-type kinetics 
derived from singular perturbation analysis of the deterministic system of equations have been used in stochastic 
simulations and the results compared to simulations for the full system [Bundschuh et al., 2003[ . There is also a 
partitioning method based on the variance of species, which led to a hybrid method coupling deterministic (small 
variance) and stochastic (large variance) by Hellander and Lotstedt ( [Hellander & Lotstedt, 2007l| ). Goutsias 
showed that fast reactions can be eliminated when probabilities of slow reactions depend at mostly linearly on 
the states of fast reactions [Goutsias, 2005[ . Utilizing the degree of advancement or extent, he followed Haseltine 
and Rawling’s approach in order to separate variables into fast and slow variables. By utilizing a Taylor series 
expansion he showed that when the slow transition rates or propensity functions depend linearly on fast extents, 
the fast reaction kinetics of a stochastic biochemical system can be approximated by the conditional mean of 
extents for the fast kinetics. Although his approach is rigorous, it has some drawbacks. In many nonlinear 
systems, slow reactions depend nonlinearly on fast reactions. For example, if a slow reaction subsystem consists 
of bimolecular reactions with two different reactants which are affected by fast reactions, his method cannot 
be applied due to the nonlinear dependence of the slow reactions on fast reactions. Peles et al applied singular 
perturbation theory and finite state projection method to obtain an approximate master equation for certain linear 
reaction systems [Peles et al, 2006[ , where clusters of fast-transitioned states have been identified. A rigorous 
mathematical framework is lacking in their derivation, which we provide herein. 

Computational methods have also been developed by many groups. Cao et al., ( [Cao et al, 2005 [ ) proposed 
a slow scale stochastic simulation algorithm (SSA). They first identified fast and slow reaction channels and 


^The notation used is defined later. 
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defined fast and slow species according as the species are changed by fast reactions or not. Assuming that the 
fast processes are stable, they approximated the propensity functions on a slow time scale, and showed that com¬ 
putation of the effective propensity functions requires first and second moments of the steady-state distribution 
of the partitioned fast reaction subsystem. They illustrated numerical results for a system with one independent 
fast species. Using a partitioning approach similar to that in [Haseltine & Rawlings, 2Q02( |Rao & Arkin, 200^ 
Cao et al, 2005 1| , Chevalier and El-Samad ( [Chevalier & El-Samad, 20091 ) derived evolution equations for both 


fast and slow reactions that led to an SSA with slow reaction trajectories in which the SSA is run until the first 
occurrence of a slow reaction. E et al. ([ |E et al., 2005| |) proposed a nested stochastic simulation algorithm with 
inner loops for the fast reactions and outer loops for the slow reactions. Strong convergence of the nested 
SSA has been proved recently by Huang and Liu [Huang & Liu, 2014| and a speed up of the the algorithm is 
proposed by using the tau-leaping method as the inner solver. More recently, Kim et al ( [Kim et al, 2014|| ) 
investigated the validity of the stochastic QSSA, where the propensity functions resulted from their determinis¬ 
tic counterpart. Under the moment closure assumption, they have shown that the stochastic QSSA is accurate 
for a two-dimensional monostable system, if the deterministic QSS solution is not very sensitive. Goutsias 
and Jenkinson ( [Goutsias & Jenkinson, 201^ ) provide an excellent review article covering recently developed 
techniques on approximating and solving master equations for complex networks. Other recent papers include 
[Cotter, 2015[|Sniith et al, 2015|| . 

As described above, many computational results and some algorithms have been reported to date, but a 
rigorous analysis of stochastic reaction networks that involve multiple time scales based on singular perturbation 
has heretofore not been done. Preliminary work on this by one of the authors [Lee & Lui, 2009[ proposed a 
reduction method based on a singular perturbation analysis for networks with two or more time scales under 
the assumption that the sub-graph of fast steps in the network is strongly connected. Our objectives here are 
to develop a rigorous analytic framework for the reduction of general stochastic networks with two widely- 
separated time scales, to prove that the generator of the reduced system is Markovian, and to illustrate the 
numerical accuracy and speedup of the reduction method for several biological models. 


2 The general formulation for reacting systems 

2.1 Deterministic evolution equations 


To set notation for the stochastic analysis, we first recall some notation used in our previous analysis of deter¬ 
ministic systems [Lee & Othmer, 200^ , hereafter referred to as I. Lurther details can be found there. 

Suppose that the reacting mixture contains the set M of m chemical species Mi that participate in a total 
of r reactions. Let un be the stoichiometric coefficient of the species in the reaction. The ua are non¬ 
negative integers that represent the normalized molar proportions or stoichiometric coefficients of the species in 
a reaction. Each reaction is written in the form 




reac a ^ . 




E prod 

ly- 


prod 

ii 


Mi 


£ = 1, ... r, 


( 1 ) 


where the sums are over reactants and products, respectively in the reaction. In this formulation, the for- 
ward and reverse reaction of a reversible pair are considered separately, as two irreversible reactions. Once 
the reactants and products for each reaction are specified, the significant entities so far as the network topol¬ 
ogy is concerned are not the species themselves, but rather the linear combinations of species that appear 
as reactants or products in the various elementary steps. These linear combinations of species are complexes 
[ Horn & Jackson, 1972[ , and we suppose that there are p of them. A species may also be a complex as is the 
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case for first-order reactions. Once the complexes are fixed, their composition is specified unambiguously, and 
we let V denote the m x p matrix whose column encodes the stoichiometric amounts of the species in the 
jth complex. 

The set of reactions gives rise to a directed graph Q as follows. Each complex is identified with a vertex Vj in 
G and a directed edge Ei is introduced into G for each reaction. The topology of G is encoded in its vertex-edge 
incidence matrix £, which is defined as follows. 


£ji — 


-hi El is incident at Vj and is directed toward it 
— 1 if El is incident at Vj and is directed away from it 

0 otherwise 


( 2 ) 


Since there are p complexes and r reactions, £ has p rows and r columns, and every column has exactly one 
-hi and one —1. Each edge carries a nonnegative weight Ri{c) given by the intrinsic rate of the corresponding 
reaction. Eor example, the following table gives four classes of first-order reactions studied in Gadgil, et al, 
([[Gadgil et al, 2005|) and two additional bimolecular reaction types. Eor either type III or VI reactions there 


Label 

Type of reaction 

Reaction 

Rate 

I 

Production from a source 

(/) Jv[i 

kf 

II 

Degradation 

A4i (/) 

kfui 

III 

Conversion 

JVlj — )' VihAi 


IV 

Catalytic production from source 

Mi 

kffn, 

V 

Bimolecular degradation 

AAj -h (j) 

kjlnjUk 

VI 

Bimolecular conversion 

Mj -h Mk 



Table 1: The four classes of first-order reactions considered in [Gadgil et al, 2005| and two types of bimolecular 
reactions. The relationship between the deterministic and stochastic rates of these reactions are discussed later. 
Rates are given in terms of number of molecules (n will be introduced later). 


may also be different types of products, e.g., A ^ B + C may represent the decomposition of a complex. 
Inclusion of such types poses no difficulties, but if such reactions are reversible we restrict the type to uni- or 
bimolecular reactions. 

If at least one reaction of type I, II, IV, or V is present the stoichiometric matrix is 

= [ I I 0 ]. (3) 

wherein the column of zeroes in the complex matrix represents the null complex 0, which by definition contains 
no time-varying species. If the system is closed and the reactions are all first order 

^ = [ I ]• 

If all species are also complexes, which can occur when there are inputs (0 ^ Mi), outputs (Mi 0), and 
first-order decay or conversion reactions (Mj Mi), the stoichiometric matrix has the form 

1/ = [ I I Z.1 I 0 ]. (4) 
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wherein vi defines the stoichiometry of the higher-order complexes. In this case the corresponding incidence 
matrix £ can be written as follows. 


£ ^ [ £i \ £2 \ £o 


(5) 


where £i represents first-order reactions, £2 represents second-order reactions, and £0 represents input and 
output steps, all of the appropriate dimensions. An alternate form of the complex and incidence matrices arises 
if the inputs or outputs are only of type I or II, for then the null complex 0 can be omitted from z/, the zbl’s 
omitted from £, and the inputs or outputs represented by a separate vector in the evolution equations given 
below. In either case the stoichiometry of the reactions and the topology of the network are easily encoded in u 
and £, respectively. 

In this notation the evolution of the composition of a reacting mixture is governed by 


dc 

- = u£R{c), c(0) = co 


( 6 ) 


where the column of 1 / gives the composition of the complex and Ri{c) is the rate of the reaction, 
or equivalently, the flow on the edge of Q. The matrix 1 ) = i/£ is called the stoichiometric matrix when 
the composition of complexes and the topology of Q are not encoded separately, as we do here [Aris, 19651 . 
One can interpret the factored form in Q as follows: the vector R gives the flows on edges due to reactions of 
the complexes, the incidence matrix maps this flow to the sum of all flows entering and leaving a given node 
(a complex), and the matrix v converts the net change in a complex to the appropriate change in the molecular 
species. 

A component is a connected subgraph Qi C that is maximal with respect to the inclusion of edges, i.e., 
if Q 2 is a connected subgraph and Qi d Q 2 ^ ^ then = ^ 2 - An isolated vertex is a component and every 

vertex is contained in one and only one component. A directed graph Q is strongly connected if for every pair of 
vertices is reachable from Vj and vice-versa. Thus a directed graph is strongly connected if and only 

if there exists a closed, directed edge sequence that contains all the edges in the graph. A strongly-connected 
component of Q (a strong component or SCC for short) is a strongly-connected subgraph of a directed graph Q 
that is maximal with respect to inclusion of edges. An isolated vertex in a directed graph is a strong component, 
and every vertex is contained in one and only one component. Strong components in the directed graph Q are 
classified into three distinct types: sources, internal strong components and absorbing strong components. A 
source is a subgraph which has outgoing edges to other strong components and has no incoming edges from 
other strong components. An internal strong component is a strong component in which edges from other 
strong components terminate and from which edges to other strong components originate. An absorbing strong 
component is a strong component from which no edges to other strong components originate. If Q has p vertices 
and q strong components then it is easily shown that the rank of £ is p{£) — p — q nChen, 19711 . 

For ideal mass-action kinetics, which we consider here, the flow on the edge, which originates at the 
vertex, depends only on the species in the complex, and the rate can be written as 


Re{c) = kejPj{c) where Pj{c) = (7) 

i=l 

for every reaction that involves the complex as the reactantj^ Thus the rate vector can be written 


R{c) = ICP{c) 


( 8 ) 


^This form also includes non-ideal mass action rate laws, but the concentrations in (|^ are then replaced by the activities of the species 
in the reactant complex, and as a result the flow on a edge may depend on all species in the system. 
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where /C is an r x p matrix with k^j > 0 if and only if the edge leaves the vertex, and k^j = 0 otherwise. 
Since each row of K has one and only one positive entry, we now denote the only positive entry in row by 

ki. 

The topology of the underlying graph Q enters into 1C as follows. Define the exit matrix of G by replacing 
all VsinE by zeroes, and changing the sign of the resulting matrix. Let K be the r x r diagonal matrix with the 
fc^’s, ^ = 1, ... r, along the diagonal. Then it is easy to see that 1C — KS^ and therefore 

j^=v£k£lP{c) (9) 

It follows from the definitions that (i) the (p, entry, p Y SKSj is nonzero (and positive) if and only if 
there is a directed edge (g,p) G G, (ii) each diagonal entry of SKE^ is minus the sum of the k's for all edges 
that leave the vertex, and (hi) the columns of EKE^ all sum to zero, and so the rank of EKE^ is < p — 1. 

If one separates the inputs, which are constants, one can write this as 

- = iy£K£jP{c) + ^, (10) 

where $ is the constant input and both u and E are modified appropriately. Herein we use the evolution equations 
in the form Q unless stated otherwise. 

One can also describe the evolution of a reacting system in terms of the number of molecules present for 
each species. Let n — (ni, n 2 ,..., rim) denote the discrete composition vector whose component rii is the 
number of molecules of species M-i present in the volume V. This is related to the composition vector c by 
n = MaV c, where A/a is Avagadro’s number, and although the rii take discrete values, they are regarded as 
continuous when large numbers are present. From Q we obtain the deterministic evolution for n as 

rhi ~ 

-^ = v£n{n) (11) 

where IZin) = MaVR{ n/N aV). In particular, for ideal mass-action kinetics 


lZt{n) 


NAVkiP^{n/NAV) 

TO / X j/j 


kj 


i=l 


i=l 


( 12 ) 

(13) 


2.2 Invariants of reaction networks 


The kinematic invariants and the kinetic invariant manifolds in a deterministic description of reactions in a 
constant-volume system are discussed in detail in HOthmer, 1979| . In general the concentration space has the 
decomposition 

nm = k[{u£f]®TZ[i^£], (14) 


where M{A) denotes the null space of A, and 72(^) denotes the range of A. The solution of ([^, which defines 
a curve in through an initial point cq, can be written 


c{t) = C{) + vE 



This shows that c(f) — cq G Tl{h'E), and the intersection of the translate of TZ{pE) by cq, which formally is a coset 
of TZ{pE) with the non-negative cone of R^, defines the reaction simplex ri(co). While this terminology 
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has a long history [Aris, 19651 , f] is a simplex in the mathematical sense only if the intersection of the coset with 
C+ is compact, which occurs if and only if there is a vector y > 0 e lOthmer, 1979| |. This is only 

guaranteed in closed systems, where the total mass is conserved and y comprises the molecular weights of the 
species, and therefore in general we should call the kinetic manifold, but we retain the standard terminology. 

First suppose that the system is closed - the case of an open system is discussed later. A vector a G Rm 
defines an invariant linear combination of concentrations if 


(a, u£R{c)) = 0, 


(15) 


and these are called kinematic invariants if a G HOthmer, 1979[ . The following important properties 

of these invariants are proven in [Lee & Othmer, 20091 . 

PI One can choose a basis for of vectors with integer entries. 

P2 If the reaction simplex ri(co) is compact, then there is a basis for M[{'u£Y'] for which all basis vectors 
have nonnegative integer entries. 


If the reactions are partitioned into subsets of fast and slow steps we can write 


i'£ — u 


£f I 


(16) 


and it follows that any invariant of the full system is simultaneously an invariant of the fast and slow systems. We 
assume throughout that the slow and fast reactions are independent, and therefore TZ{'u£) — TZ{v £^) U TZ{'u£^). 
Furthermore, one has that 


C (17) 

V[(j/^)^] CA^[(z.f(/’*))^]. (18) 


that the ranges of the slow and fast subsets are no larger than that of the full system, and the corresponding null 
spaces are no smaller, but the properties PI, P2 of the full system do not necessarily carry over to the subsystems. 
However, it follows from (18) that one can define a map ^ Rm-r/ (where Vf - AimR[v£f ]) for the 

fast subsystem that represents a vector in M[{'u£f)'^] in terms of intrinsic coordinates on M[{'u£f)'^]. The 
associated matrix has rows given by basis vectors with integer components of It follows that 

the reaction simplex for the fast subsystem is given by 


f]j(co) = {c : c G Co + TZ[v£^]} n R^ = {c : c = 'P-^cq = c G H R^ 


Here c represents a conserved quantity for the fast subsystem, but it may vary as slow reactions occur. 

If the system is open then one can regard the effect of inputs and outputs as moving the dynamics between 
simplexes of fixed mass, as dictated by the evolution equations given at Q. This point of view will be useful in 
the stochastic formulation. 


3 The stochastic formulation 

3.1 The master equation 

In a stochastic description the number of molecules of a species is too small to be treated as a continuous variable 
- they are random variables. We define N{t) = (ni(f), n 2 (f), ..., where ni{t) is as before, but now 
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N(t) is a random vector. Under the assumption that the process is Markovian, which is appropriate if all the 
relevant species are taken into account, the evolution of N (t) is governed by a continuous-time Markov process 
with discrete states, and we denote the probability that {N{t) = n} by P(n, t). The governing equation for the 
evolution of P(n, t) is called the chemical master equation, and is given as 

^P(n, i) = Uiin - vS(i)) ■ P{n - Uein) ■ P{n, t) (19) 

e e 

where the column of 6, denotes reaction and the stochastic rates have the form 


TZi = (20) 

Here C£ is the probability per unit time that the molecular species in the complex reacts, j (£) denotes the 
reactant complex for the reaction, and (n) is the number of independent combinations of the molecular 
components in this complex [ Gadgil et al., 2005] |. One sees in ( [T^ how the reaction graph Q determines the 
state transition graph Qg that defines the steps in the master equation. 

In the stochastic analysis all reactions are assumed to follow mass-action kinetics, and thus q = ki, and the 
combinatorial coefficient is given by[^ 

nC";V 


If = 1 the stochastic rate reduces to (13) but for a bimolecular reaction the rate of a step in the stochastic 
framework is always smaller than in the deterministic framework. 

The master equation ( [T^ can be expressed more explicitly for unimolecular and bimolecular mass action 
kinetic mechanisms. The general form for the uni- and bimolecular molecular reactions given in Table[^is 

dP{n^ t) 


dt 


= E 

i=l 


+ E 

k = l 


Kt,{Sr^ - l)P{n,t) + ^ - 1) + K^\S-^ - 1) + (5+^ - 1) 

J = 1 ^ 

- l) + KfiiS+^St - l)y,ynjPin,t)) . 


where Sf is the shift operator that increases the component of n by an integer amount k, and Sf {niP{n^ t)) = 
Sfrii • P{Sfn^ t). In this expression the matrices and correspond to the six 

types of reactions given in Table 

A compact representation of ( [T^ is obtained by defining an ordering of the accessible states, of which we first 
assume there finitely many (n^) in number and which lie in the non-negative cone C+ C The evolution 

of the vector p{t) of the probabilities of these states is governed by the matrix Kolmogorov equation 

dp{t) 


dt 


= Kp{t), 


( 22 ) 


where K is the matrix of transition rates between the states, the entries of which are defined as follows. Let 
iP = {n \, 77.2,, ^m) — (^i 5 ^ 2 ’ • • • ’ denote the and states of the system and denote the 

and entries of the vector p as = P{rd^ t) and pj = P{n ^, f), resp. Then the (i, entry of K is given by 




if n'^ = for some ^ = 1,..., r 

0 otherwise 


(23) 


^This formulation applies only to ideal solutions - in nonideal solutions the number of molecules must be replaced by an appropriate 
measure of its activity in the solution. In particular, this involves a suitable description of diffusion when the solution is not ideal 


lOthmer, 1976|[Schnell & Turner, 2004) . 


The case in which there are infinitely many states will be discussed in Section 3.2 


8 














and the following pseudo-code shows how to generate K and the corresponding state transition graph Qs foi* 
any reacting system. By starting with a given initial state vector, the algorithm checks all the reactions and 
adds resulting new state vectors into state space while updating the transition matrix, and then repeats the same 
procedure on the newly added state vector until no new state vector can be added. 


Data: initial state vector Vi and reactions 
Result: transition matrix 

Initialization: set current state index = 1, set accessible state vector space {!/} to be {Vi}, set accessible state 
vector space size = 1, set transition rate matrix K to be {0}; 
while current state index < ris do 

set current reaction index to one: = 1; 

while current reaction < r do 

check if current reaction Rc reacts from current state Vc ; 
if True then 

get current state index: source ^ Cs\ 
get target state Vt = Vc ^ ’ 

check if target state Vt is already in {V}; 
if True then 

get the index i of the state in {y} that is equal to Vt; 
check if Ki gQy^Y^QQ — 0, 

if True then 
I update 
end 

end 

else 

add Vt to {V}: K, ^ Vt; 

increase accessible state space size by one: Ug ^ -h 1; 
update Kji^^gQy^j^QQ, 

end 

end 

increase current reaction index by one: ^ -h 1; 

end 

increase current state index by one: Cg ^ Cg + 1; 

end 

update diagonal entries Kjj = — i^j Kij. 

Algorithm 1: An algorithm for generating the transition matrix K for a chemical reaction network 
Off-diagonal elements in the row of K are transition rates at which source states reach state i in one step, 
whereas off-diagonal elements in the column represent the rates at which the state reaches its target states 
in one reaction. Since only uni- and bimolecular reactions are realistic, the states that reach the current state in 
one step (the sources) differ from the current state by at most two molecules, but those that can be reached in 
one step (the targets) may involve more than two. Under certain orderings K has a narrow bandwidth when the 
number of states Ug is small, but expands with increasing Ug. The number rig grows combinatorially with the 
number of molecules - when there are m species and a total of Nq molecules 

_ /TVo + ^ - A 
^ \ m — 1 J 


If, for example, there are 4 species and 50 molecules, the number of states is 23,426. In any case, if all reactive 
states are accounted for the matrix K is the generator of a Markov chain [Norris, 1998|| , i.e. ^ij — 0 
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each j, Kij > 0 for each j ^ i, and the vector = (1,1,, 1)^ of length equal to rig is in 
The formal solution of ( [22| ) is 

p{t) = e^V(O), 

but there are well-known difficulties in computing the exponential if the state space is large, and various approx¬ 
imate methods have been used [Kazeev et aL, 2014] |Menz et al, 2012[ |Deuflhard et al, 2d08l |. An alternate 
approach is to do direct stochastic simulations, using the Gillespie algorithm or one of its many variants. How¬ 
ever the computational time can be extremely long if there are many species or the system is spatially distributed 
iHu et ai, 20131 , and thus it is advantageous to reduce the system if possible. 

Reduction of deterministic systems is frequently done by taking advantage of the presence of multiple 
time scales in the evolution, and reducing the number of variables by invoking the QSSH, as was done in 
|Lee & Othmer, 200^ . In the stochastic system one can determine whether a reaction is fast or not according to 
the magnitude of the transition rate of a reaction, the so-called propensity function in the chemical literature. A 
larger transition rate implies that the corresponding step occurs more frequently, but since the transition rate of 
a reaction depends on the numbers of reactant molecules as well as the rate constant of the reaction, a reaction 
that is fast in the interior of C+ may be slow near the boundary of the cone. The following example shows that 
a reaction with a large rate constant may have to be considered as a slow reaction. 


Example 1 Consider the following reaction network 

K 

B 

k-i 

Ifki = 0.01, k-i = 0.1, and initially 72^(0) = 100, ^^(0) = 1, then the rate of A B is kiriAf^) = 1, and 
from B ^ A is k-insiO) = 0.1. Thus even though ki < k-i, initially A ^ B is fast compared with B ^ A. 
Of course when A is small kinA{t) < k-insit), and which is fast and which is slow is interchanged. 


This issue must be addressed on a case-by-case basis, but here we simply assume that all steps can be identified 
as occurring on either an O(^) scale or on an 0(1) scale a priori. Since we consider a small system with 
small numbers of molecular species involved in at most second-order reactions, the separation assumption can 
be justified by stipulating that 

kf « of) and N^ks « 0(1), 

where A^o is the maximum number of molecules of reactants involved in slow reactions and kf and kg are 
characteristic rate constants of fast and slow reactions, respectively. Thus we can rewrite the matrix K as 

K = -Rf + K^, 
e 


where and are the fast and slow transition matrices whose entries are the transition rates of fast and slow 


reactions, respectively. Then on the 0(1) timescale (22) can be written 


dp{t) 

dt 


{ff + K^)p{t). 


(24) 


As before, conservation of probability must be satisfied, but now it must hold separately on the slow and fast 
scales, i.e., G and A/^((iT^)^). Thus both the transition rate matrices and for slow and 

fast reactions are generators of Markov chains. As we see later, the sub-graphs qI and Ql of Qg associated 
with the fast and slow reactions may have several disjoint connected components when considered as undirected 
graphs. These need not be strong components of Qg. 
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Equation (24) is written on a slow time, but on the fast timescale r = f/e it reads 


^ = + eK^)pit) (25) 

and from this one sees that the slow reactions act as a perturbation of the fast reactions on this scale. It follows 
that the perturbation can only have an effect on the dynamics represented by the eigenvector(s) corresponding 
to the zero eigenvalue of K^, and this will be exploited later. 


3.2 Conditions under which the state space is bounded 


A number of technical issues have to be discussed before proceeding further. Firstly, it is easy to see that the 
positive cone of is invariant, which therefore preserves non-negativity of the probabilities. Furthermore, 
if the reaction simplex is compact there are only finitely many states in the system, and the generator of the 
Markov chain is a finite matrix. In that case the invariant distribution is well-defined and, given certain additional 
properties discussed later, it is unique. When the state space has infinitely many points the generator is an infinite 
matrix, and in fact, is unbounded as an operator on /^xj^because the entries of K depend on n. In this case little 
is known in general about its spectral properties and invariant measure(s). This is in contrast with random walks 
on homogeneous lattices, where the generator is bounded even if the state space is infinite. Moreover, the fact 
that the corresponding deterministic description has a compact invariant set is of no import in the stochastic 
description, since large deviations from the mean dynamics are possible. However, as the following example 
shows, the probabilities of very large numbers may in general be very small under suitable conditions. 


Example 2 Consider the simple process (j) A 0, where as usual, (j) represents a source and sink, and 
let Pn{t) be the probability of having n molecules of A in the system at time t. Then one can show that the 
generating function 


Gis,t) = yy^s^pnit) 


n=0 


satisfies 


dG , , ^.dG , , 

— + k2{s-l)—-h{s-l)G. 


The solution of this for the initial condition pi{0) = l,p^(0) = 0 otherwise, is 
G { sA ) = (^1 + (5 - “ 1)(1 “ 


and from this one finds, by expanding this as a series in s, that 

n—l 

' I _ g-fet 


i (I) (‘ 

Therefore the stationary distribution is 


^ (^ki{l — e -h k 2 ne • exp 


1^1 


—k2t^ 


(26) 


1- / ^ 1 

hm p„(f) = — — exp 

t^oo n\ \k 2 J 


h 

'k2 


which is a Poisson distribution with parameter ki/k 2 . Thus Pn{t) is non-zero for arbitrarily large n in both the 
transient and stationary distributions, but it decays rapidly with n. For example, ifki/k 2 ^ 0(1) (^kid n = 25, 
Pn ^ 0(10“^^) in the stationary distribution. Even if the stationary mean k\/k 2 ^ 0(10), Pn < for 

n 50 (one must always choose n greater than the mean in order that pk < Pnfar k > n). 

vector space whose elements are infinite sequences of real numbers. 
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While this is only a heuristic justification for truncating the state space to a box of the form [0, Mi ], we 
shall do this in the remainder of the paper. One could formalize the truncation by modifying the dynamics so as 
to make all states that lie in the positive cone and outside the box transient, and then restricting the initial data 
to lie in the box. Further work on the validity of this truncation is needed, but in any case the following results 
apply rigorously to closed systems and those systems whose inputs are terminated when the total number of the 
molecules being input reaches a pre-determined level, as in the following example. 

Example 3 Consider a triangular system in which there is a input of one species, as shown below. 

k2 k4 ke 

If the initial State is (72^(0), 77 . 5 ( 0 ), nc'(O)) = (0,0,0) this defines the initial simplex by nA{0)-\-n B{0)-\-nc{0)> 
Whenever a molecule of A is added the simplex shifts upward by one unit in the first octant. To insure that the 
simplex remains bounded we simply turn off the input f when the number of species A achieves a certain 
threshold, say, Na> At this moment t — C, we also observe the molecular number of species B and C, which are 
nB{C)^Tic{C), and the system evolves on this reaction simplex thereafter, since the system is closed thereafter. 

In general we shall assume that the inputs and outputs occur in the slow dynamics, and thus the fast reactions 
occur on a fixed finite simplex. The inputs and outputs can move the system ‘up’ and ‘down’ in the positive 
cone, but we assume that this occurs on the slow time scale and if the total inputs are controlled the state space 
will remain bounded. 


3.3 The role of invariants 


Any conservation condition that applies to a system defined in terms of concentrations applies when the system is 
defined in terms of numbers, and therefore the invariants that characterize Q apply in the stochastic framework 
as well. In particular, the definition of a kinetic manifold Tt f for the fast subsystem, in which the slow variable c is 
invariant, carries over with only minor modification. Since the state space is discrete in a stochastic description, 
only the set of integer points in a reaction simplex Q are relevant, and we call the subset of integer points in Q 
the discrete reaction simplex and denote it by jC. 

For any c G ^^(co), we can write c = cq + oSp > 0, where p is the extent, and therefore 

n = JVaVc = JVaVcq +AfAVoSp G Q{J\fAVco) =: ^^(77.0). 

Since n lies in Z^, which is the closure of the set Z+ of 777.-dimensional vectors with positive integer entries, it 
follows that the set of all integer points in the reaction simplex is given 

£(no) = f^(no) n Z^, 

and we call this the (full) discrete reaction simplex through 77 . 0 . Analytically, a discrete simplex jC{no) is defined 
as a coset of IZfuS), and it can be generated numerically using Algorithm To illustrate this, consider the 
following reaction network for an open system. 


Example 4 Consider the network 

(j) ^ A^ B ^ (j). 

where the input and output reactions are slow steps. One finds that 


TZ{o£) n Z^ = span 
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which is Z^. If we start with the initial state (1,0), and stop the input reaction once one molecule of A is added 
from the source, as described in Section ^?^ we can generate the full discrete reaction simplex as the set of all 
points denoted Figure^a). 

Since the inputs and outputs occur on the slow time scale, these reactions are denoted by red arrows in 
Figure 3.1 (b), and the remaining reactions are in fast components (green arrows). There are two distinct 
fast (absorbing) strong components, as shown in Figure^b), defined as C{ = {(1,0), (0,1)} and — 
{(2,0), (1,1), (0,2)}. Since 7l(uE^) = span{[—l = span{[l 1]}, we can identify each 

fast simplex uniquely by 


h = A^n — [1 1] 


ni 

n2 


[ni +^ 2 ], 


where ni{n 2 ) denotes the number of molecules of species A (B). The rows of the matrix A^ comprise a basis 
for and thus C{ is defined by h — 1, and C 2 by h — 2. In this simple example, we can interpret h 

as the total mass, which is conserved by the fast reactions. 


B B 




Figure 1: (a) The full discrete reaction simplex; (b) two fast discrete reaction simplexes distinguished by yellow 
and blue colors. 


To define fast simplexes in general we consider only fast reactions, and first identify all distinct absorbing 
strong components in qI . Then, working backwards, we identify the sources and internal strong components 
that feed into an absorbing component. This identifies a sub-graph of Qs comprised of an absorbing strong 
component and its ‘feeders’, and leads to the following definition of a fast simplex. 

Definition 5 A fast discrete reaction simplex is the set of vertices in an absorbing strong component and its 
predecessor sources and internal strong components. 

Remark 6 

1. A fast simplex is uniquely determined by its absorbing strong component. Each fast simplex is character¬ 
ized by the invariants of the fast dynamics. These are 

h = A^n (27) 

where the rows of the matrix A^ comprise a basis for )^]. 
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2. The invariant distribution of the fast dynamics on a given fast simplex is non-zero only for the states in the 
corresponding absorbing component. 

3. Both source components and internal strong components may belong to more than one fast simplex, as is 
exemplified by a reaction network with a binary tree structure. In such a network the root belongs to all 
fast simplexes and internal nodes at the level belong to simplexes, where n is the depth of the 
tree. 

4. If there are no slow reactions in the network there is only one fast simplex, and this case is not of interest 
here. 

5. In general the disjoint union of all fast simplexes is the full fast simplex, i.e. 

= Vtf (no) n Z+. 

Although the following reaction network due to Wilhelm ([ [Wilhelm, 2009| ) involves unrealistic trimolecular 
steps, we use it to illustrate how to define the fast simplexes when the fast subgraph qI has multiple absorbing 
components whose predecessors are not disjoint. 

Example 7 


S + 2X 2.x 

3X 2X + P 


The species are S, X, P, and the complexes: S+2X, 3X, 2X-\-P, X, P are labeled 1-5 in this order. Then the 
matrices v, £ and v£ are as follows. 







-1 

0 

0 





’ 1 

0 

0 

0 

0 ' 


1 

-1 

0 


’ -1 

0 

0 

V — 

2 

3 

2 

1 

0 


0 

1 

0 


1 

-1 

-1 


_ 0 

0 

1 

0 

1 


0 

0 

-1 


0 

1 

1 






0 

0 

1 





It follows that p{£) = 3^p{o£) = 2, and the deficiency, which is defined as 5 = p{£) — p{o£) = 1. 
Clearly S -\- X -\- P is a constant and equivalently, M[{o£Y'] — span ((1,1,1)^). If we set the initial state as 
(n5'(0), nx(0), np(0)) = (1, 2,0), then the possible transitions in the system are as shown in Figure^a). 

Assume that only the second reaction is slow and the others are fast. Considering only the fast transitions 
indicated by green directed edges, it is clear that the resulting qI has the two absorbing strong components 
(1, 0, 2) and (0,0, 3) because there are two possible exit steps from the node (1, 2,0), one with rate ki, the other 
with rate k^. The two fast simplexes are 

4 ={( 1 , 2 , 0 ),( 1 , 1 , 1 ),( 1 , 0 , 2 )} 

4 ={( 1 , 2 , 0 ), ( 0 , 3 , 0 ), ( 0 , 2 , 1 ), ( 0 , 1 , 2 ), ( 0 , 0 , 3 )} 

as shown in Figure^b) & (c). Figuratively, and also in the language of graph theory, the simplexes are branches 
of the tree that defines the full simplex in Figure^a). jC{ lies on the coset oflZfvS ^^) = span{ [0 — 1 1]^}, 
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Figure 2: Examplej^ (a) The full discrete reaction simplex(black circles); (b) & (c) the fast simplexes. 

thus — span{[l 0 0],[0 1 1]}, and C{ is identified by hi = (1, 2). lies on the coset of 

7l{iy£^^) = span{[—l 1 0]^, [0 — 1 1]^}, thus = span{[l 1 1]}, and jC^ identified 

by 02 = 3. 

In general a component of the reaction graph may have more than one absorbing component and more than one 
fast simplex, as in the foregoing example. To determine them in general one must first identify the number of 
strong components and their type in every component graph Ga ^ G • Since no vertex in a source is reachable 
from any vertex outside its component and no vertex in a sink is reachable from a vertex in any other sink, the 
relationship of reachability defines a partial order on the strong components of Ga- This in turn leads to the 

o o 

acyclic skeleton Ga of which is defined as follows. Associate a vertex Vj with each strong component of 

o o 

Ga, and introduce a directed edge from Vi to Vj if and only if one (and hence every) vertex in Vaj is reachable 

o 

from Vai. Ga is connected since Ga is connected, but it is acyclic; in fact, it is a directed tree. One of two 

o 

cases obtains: either Ga consists of a single vertex and no edges, which occurs when Ga consists of one strong 
component, or it has at least one vertex of in-degree zero and at least one vertex of out-degree zero. One can 

o 

relabel the strong components if necessary, and write the adjacency matrix of Ga the form 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

X 

X 

0 

0 

0 

0 

X 

X 

X 

0 

0 

0 

X 

X 

X 

X 

0 

0 

X 

X 

X 

X 

0 

0 


where the x’s represent blocks that may be non-zero. The diagonal blocks are square matrices of dimensions 
equal to the number of sources, the number of internal strong components, and the number of sinks, respectively. 
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The vertices corresponding to internal strong components can always be ordered so that the central block is lower 
triangular, since the strong components are maximal with respect to inclusion of edges. The number of sinks or 
absorbing strong components is the number of zero columns of A. 


4 The reduction of a master equation with two time scales 

4.1 The splitting of the evolution equations 

In this section we show how to obtain the lowest-order approximation to the slow dynamics for a system that is 
described by 

| = + (28) 

on the 0(1) timescale and by 

+ eK^)p (29) 

dr 

on the fast time scale r ^tje. 

By a suitable ordering of the states we can write the fast transition matrix as a block diagonal matrix 




K 




Ki 


Rf 


R 


(30) 


wherein the number of blocks is equal to the number of fast components in the state transition graph Ql of the 
fast reactions. We show later in Theorem[TT]that each block has zero as a semisimple eigenvalue and we analyze 
the structure of the blocks in detail, but here we first define the index of a matrix [ Campbell & Meyer, 1991[ . 


Definition 8 Let A : ^ be a linear transformation. The index of A is the smallest non-negative integer 

k such that 7^(A^+l) = 7^(A^). 


An equivalent definition is that the index of A is the largest Jordan block corresponding to the zero eigenvalue 
of A. It follows from the general theory in [ [Campbell & Meyer, 1991| j that a matrix with a semisimple zero 
eigenvalue has index 1, and as a result, the following properties hold. 

(i) Rn^JV(A)e7Z(A) g 

(ii) 7Z(A^) = 7Z(A) 


Since we assume that the inputs and outputs only occur on the slow time scale, the fast dynamics occur on a 
fixed compact simplex in C+ and R-^ is a bounded operator that generates a Markov chain, since 1 G J\f{R^ )^. 

^This is the direct sum, but is generally not the orthogonal direct sum. 
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We define dim{J\f{K^)) = rif, and of course it follows that both and have a basis of rif 

linearly-independent eigenvectors. Let II be the rig x rif matrix whose columns are eigenvectors corresponding 
to zero eigenvalues of K^, and let L be the nj x rig matrix whose rows are eigenvectors corresponding to zero 
eigenvalues of {K ^)^. One can choose these to form a biorthogonal set and therefore Lfl = Furthermore, 
riL = Po, where Pq is the eigenprojection corresponding to the zero eigenvalue of , and thus P(Po) = 
M{Kf) and = P(Po) © P(/ - Pq). 

By property (i) above = Af{K^) © TZ{K ^) and to show that TZ{I — Pq) = TZ{K ^) we observe that the 
adjoint P^ leads to the decomposition R^^ = P(Po^) © P(/ — P^). We have that 71{Pq) = is 

orthogonal to P(/ — Pq), and from the basic properties of projections and their adjoints |Kato, 1966| it follows 
thatP(/-Po) =n{Kf). 

Since R^, = Af{Kf) © we can write 


p = Pop + (/ - Po)p = lip + Tp. 


(31) 


where p — Lp and F is an x (n^ — rif) matrix whose columns are basis vectors of P{Kf). Then the matrix 
[n I r] is an ng x Ug matrix whose columns are basis vectors of Rn^- Using (|3^ in (|28|) leads to 




(32) 


To obtain the evolution equations for p and p separately, we define the {rig — rif) x rig matrix V = 
[0 I Ins-nf][P- I which has the property that V[Il | T] = [0 | Im-uf]- In other words, VH — 0, VT = 
and by multiplying (32) by L and V, respectively, we obtain 


^ = LiF^(np + rp), 

^ = Vi^Rf + K^){Up + Fp) = ^VRf Fp + VK^iUp + Fp). 

where we used the fact that LF = 0 because V((iF-^)^)) -L P{K^). 

On the fast time scale t = tjevi/e obtain in the limit e ^ 0 that 

^ = VK^Tp. 
dr 


(33) 

(34) 


(35) 

(36) 


On this time scale the slow variable p remains constant while the fast variable p evolves according to ( [36| ). 
On the slow time scale 

dp 


= LiF*(np + Fp), 

= {VRf + eVK^){U.p + Tp). 


(37) 

(38) 


and the limit e ^ 0 leads to the condition that 


VK^Vp^Q. 


Lemma 9 


only if Fp = 0. 


VK^ Fp = 0 
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Proof 1 By definition Tp G TZ{K^), and since — TZ{K^) by property (ii) above, it follows that 

K^Tp = Tqfor some q. Since VT = VK^Tp = 0 only ifq — 0, which implies that Vp — 0. 

Since Tp = 0, it follows from the equation ( [T7| ) that 

^ = LK^Up = Kp, (39) 

dt 

which shows how the slow evolution between fast simplexes depends on components of the fast evolution 
through L and II. The components of p represent the probabilities of aggregated states, and since Tp = 0 
it follows from <0 that p = Tip, from which one obtains the probabilities of the individual states. 

The following example illustrates the structure underlying the general reduction. 

Example 10 Consider a closed triangular reaction system with all reversible reactions, in which reactions 1 
and 2 are fast, as shown in Figure^ The total number of molecules is conserved in a closed system of first order 
reactions of this type, and the transitions in state space are as illustrated in Figure^when the total number of 
molecules is two . 



Figure 3: Triangular reaction: red and green arrows denote slow and fast reactions, respectively, (a) The reaction 
network; (b) The state transition diagram for a total of two molecules. 


The matrix is obtained by setting the rates of all slow steps to zero. The resulting graph has three 

connected components, (2,0,0)^ ^(1,1,0)^ ^(0,2,0), (1, 0,1)^_^(0,1,1) and {0^ 0^2), and this leads to 

the representation 


K = - 
e 


' k( 

0 

0 


1 

^1,2 

^1,3 

0 

-^2 

0 

+ 

1—1 


^2,3 

0 

0 

Ki . 


. ^3,1 

^3,2 

. 


(40) 


where K- fi = 1,2,3 the rrii x rrii matrix of transition rates within each fast component, and Kf ■ is an 
rrii X ^7 matrix of slow transition rates between fast components. In this example, the Kl are 


k( = 


-2ki 

2ki 

0 


k2 

-{ki + k2) 

h 


0 

2/c2 

-2k2 


1- 

-fci 

k2 


ki 

-k2 


Ki = 0 , 


Note that the vector 1 of the appropriate dimension is a left eigenvector of each KI, and therefore probability 
is conserved on each fast component. In view of the block structure of in {40) zero is a semisimple eigenvalue 
ofKf. 
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The sub-matrices of are 


—2kQ 0 

0 ~(^3 + ^e) 


0 


0 


0 

0 

-2^3 


— 


-{k^ + ^5 + k^) 

0 


0 

— {ks + ^4 + k^) 


,KI = -{2k^+2k5) 


and 


Kii = 


2kQ ks 0 
0 fee 2 k^ 


-1 


0 

, K12- 

^4 

k^ 

- 

0 

k/^ 


1 - 1 

0 

0 

0 

, KU = 

0 

0 

0 

1 1 

T 

KI2 = 

ke ks 


2/^5 2/^4 


Kh = 


The matrices L and E are given by the following. 



Li 

0 

0 


ni 

0 

0 

L = 

0 

L2 

0 

, n = 

0 

n 2 

0 


0 

0 

Ls _ 


0 

0 

ns 


where Li — [11 1], L 2 — [1 1], L 3 = 1 are basis vectors ofJ\f{Kf The E^'^ are the multinomial invariant 
distributions on the fast components l[Gadgil et al., 2005]/ , and are found to be 


Ei = 


-1 

to to 

2 k\k 2 


T 

Eo — 

k 2 

ki 

iki + k 2 V 

{ki + k 2 Y 

(fci + k 2 Y . 

5 

^^2 — 

ki + k 2 

ki + k 2 


E 3 = 1 . 


O’s are zero row or column vectors with the appropriate dimensions. 

As shown above, the dynamics on the slow time scale are governed by ( |j9| ), and in this example the transition 
rate matrix for the slow system is given by 



■ LiKfUi L 1 KI 2 U 2 LiKf ^Us ■ 


kl 

Ls 

^ 1,2 

Z.S 

^1,3 

K = LK^U = 

L 2 Kl^Ui L 2 KIU 2 L2KI:,U3 


Z.S 

^ 2,1 

Z.S 

^2 

Z.S 

^2,3 


LsKl^Ui L3Kl2^2 L 3 KIU 3 


ILs 

L ^3,1 

H ,2 

k . 


where the off-diagonal transition rates kf ■ are given by 


K 2 = LrKlffl2. 


^1,3 


— ^3,1 


= 0 , 


kh = L 2 K, 




^2,3 


= L2i^|,3n3, 


fc|,2 = ^3i^|,2n2. 


The diagonal elements are the negatives of the corresponding column sums and therefore K is the generator of 
a Markov chain. We will show later - in Theorem\T3\- that the transition rate matrix for the reduced system is 
alway a Markov chain generator. Recall as discussed in Subsection\4.3\ we have 


-1 1 

1 -1 


1 1 0 

h — A^n = 

111 + n 2 

_ 0 0 _ 


0 0 1 


m 
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Thus hi — (2, 0), h 2 — (1,1), and hs = (0, 2). The evolution on the slow time scale is as shown schematically 
below. 

Ls Ls 

^2 1 ^3 2 

( 2 , 0 )^i—►(!,!)( 0 , 2 ) 

Ls Ls 

^ 1,2 ^ 2,3 

The matrices for the general case in which the total number of molecules is A^o are given in the Appendix. 


4.2 The structure of the generator of the slow dynamics 


We recall the standing assumption that there are no inputs or outputs that occur on the fast time scale, and 
therefore on the fast time scale the reactions are confined to fast simplexes as defined earlier. The slow dynamics 
moves the state between fast simplexes 12/ on a given simplex 12, and/or between simplexes on the slow time 
scale. Because the transition matrix for the fast dynamics has the block diagonal form given in ( [30| ), both L and 
n have the same block structure and it suffices to determine the null spaces for a fixed block. 

The graph can be decomposed into sources, internal strong components and absorbing strong compo¬ 
nents. If all fast reactions on fast simplexes are reversible, then all fast components are strongly connected. If 
a state or vertex is not connected to any other states, we call it an isolated state, and it is then a (absorbing) 
strong component itself. Note that in the triangular reaction network as in Example [T0| each fast component is 
an absorbing strong component. 

We let Si be the set of all states (nodes) in the graph of the simplex and denote the set of the states in the 
sources, internal strong components and absorbing strong components by and respectively. The 

set of all the states, S can be represented by disjoint unions 


We assume that and = m|^. Thus rrii = mf^ + + m| 


^so 

H • 


By fol 


kI in (|30|) as an upper triangular block matrix 


owing the general analysis of reaction networks developed in [Othmer, 1979| | we can write each block 


-p^f^ah^in j^f,ab,so 

where each block is the transition matrix between states where a = ab^in^ so and is the 

transition matrix from the states in into states in Sf, where a, /3 = m, so or ab. 

If the simplex includes pi isolated absorbing states and qi absorbing strong components with at least two 

states, then the matrix can be written as 



7 

T^f,ab 

i^qi -I 

where 0^. is a x pi zero matrix that reflects the transition rates between the pi isolated absorbing states and 
is an x transition matrix describing the transitions between all states in the absorbing strong 
components with at least two states. 
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The following result defines an important spectral property of each which determines the structure of 


its left and right eigenvectors Li and as shown in Subsection 4.3 


Theorem 11 Each matrix has a semisimple zero eigenvalue and each and i = 1, 

nonsingular. 


Proof 2 The proo f o f this follows from the general result given in HOthmer, 1979]/. This theorem implies that the 
zero eigenvalue of K- is semisimple. 


4.3 The invariant distributions of the fast dynamics 

Next we consider the invariant distributions of the fast dynamics, which give the basis vectors in II. These are 
vectors tt > 0 such that — Q and Hj — 1. Since every absorbing strong component with at least two 
states has a unique steady-state probability tt^^ , j — 1, • • • , is a basis for M since — 0, 

and diva M Notice that each absorbing strong component with only one state, which corresponds to 
each diagonal entry of the block 0^. in has steady-state probability 1. Using these facts we define 




^Pi 


7T, 


ah 

z,l 


.ah 


where Ip. is an Pi X pi unit matrix and irfj is an x 1 vector. Since any states in sources and internal strong 
components have zero probability at the steady-state, we define the matrices 



n* 


ni 


O 

O 

o 




Osi Osi ■ * * Osi 




where and O^. are null matrices representing the steady-state probability of all states in internal strong 
components and sources, respectively. Since K^Ui = 0 it follows that K^H = 0, and therefore we have an 
orthogonal basis for 

To construct a basis of )^), recall that the multiplicity of an eigenvalue is the same for a matrix and 

its adjoint, and therefore has exactly one zero eigenvalue and so dim(A/^(iT/^^^)^) = 1 for each j. The 

1 X vector = [1,..., 1] is an eigenvector corresponding to the zero eigenvalue of N, and 
therefore is a basis of for each j. Define 


^Pi 
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and note that LiKf’°'^ = 0. By defining Lj = Lj Ai Bi 
Bi = one can see that 


, where Ai = and 


LiK! = 


Lj, Bi 


f,ab j^f,ab,in j^f,ab,so 


Kr Ki 

K: 


f,in 


K, 


f,in,so 


K: 


f,so 


= 0 . 


Thus, 


where 


LRf = 0 , 


L = 


(42) 


One can see that all row vectors of Li consist of a basis of N{{KfY') and it follows from the block diagonal 
structure of L that all row vectors of L are basis vectors of One can also show that the sets {Lj} 

and {Ilj} are biorthogonal. 

Remark 12 When all components are strongly connected, e.g., all fast reactions are reversible, the matrices L 
and n reduce to 



Li 


Hi ... 

L = 


, and n = 



Li 


n; 


where Li is al x rrii vector [1,1 • • • 1] and is anrriiX 1 steady-state vector of a strong component Li. Here 
rrii is the number of states in Li. 


The following theorem shows that the slow evolution is always Markovian, even when a fast simplex has 
multiple components, each of which may have sources and/or internal components. 

Theorem UK — LK^Ii is the generator of a Markov chain. 


Proof 3 In order to prove this, we only need to show that the sum of each column of K is zero and the off- 
diagonal entries are nonnegative. Without loss of generality and to avoid complicated subindices, we consider 
the first column of K. The first column block of K, denoted as is given by 

{k)^ = [LrKiBr L2iT|,ini ••• 


where n is the number of fast components, block G Llabixabifi — 1 ,... abi is the number of 

absorbing strong components of the i-thfast component. The sum of all the rows of Ki is 


i ^ ki ' i ki i ^ ki 

— llxDi-fi”**,!!!! = ( y^ = OixDiHi = Oixafci, 
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where denotes the ki-th row of matrix A and Di is the number of nodes in the i-th fast component. The 

third equality follows from the fact that the column sum of Li is always one, because each Li defines an invariant 
probability and the sum must be one. The last equality follows from the fact that the column sum of is zero. 
Since non-positive entries appear only along the diagonal of Kf the off-diagonal entries of K are nonnegative. 


To summarize, we now have the approximate equation for the probability distribution on the slow time scale 


dp 

dt 


LK^Iip, 


(43) 


where the matrices L and E are given above. Clearly the utility of the reduction depends heavily on whether 
the invariant distributions of the fast dynamics are easily computed. This is possible for a network of first-order 
reactions, and for closed systems the distribution is multinomial if the initial distribution is multinomial, while 
for open first-order systems the distribution is a product Poisson [Gadgil et al, 2005| . The general case of first- 
order reactions is completely solved in [Jahnke & Huisinga, 2007] |, where the authors show that the solution 
can be represented as the convolution of multinomial and product Poisson distributions with time-dependent 
parameters that evolve according to the standard equations. 

Analytical results for bimolecular reactions are more limited. However it is known that the invariant distribu¬ 
tions of a class of Markov chains that arise in queuing theory have a product form [Boucherie & Dijk, 199 1| , and 
it has been shown that this also holds for a zero-deficiency reaction network (which means that TZ{£) = 

(/)) under the assumption of ideal mass-action kinetics [Anderson & Kurtz, 201 1[ Melykuti et al., 20f4| . When 
this does not apply the stationary distributions frequently involve hypergeometric functions [McQuarrie, 1967 1. 


4.4 A detailed example of the reduction 

To illustrate the reduction method on general reaction networks, we consider the three reaction networks shown 
in Figure]^ The resulting states beginning from one molecule of A and two of F are shown in Tablewhere the 
states are in the order of molecular species A, S, • • • , G. The state space of network 3 is a subset of that of net¬ 
works 1 and 2 since reactions with rates and k^ are missing, as seen in see Table(right). The corresponding 
full state networks as shown in Figure]^ are generated using MATLAB. 
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Table 2: Accessible states of networks. Left: Network 1 & 2; Right: Network 3 

By switching off the slow reactions in the graphs of the full state transition diagram shown in Figure]^ we 
obtain three fast components for each network as shown in Figure]^ (top), (center) and (bottom) respectively. 
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Figure 4: (a) Network 1. (b) Network 2. (c) Network 3. Red and green arrows denote slow and fast reactions, respectively. 


The fast components on the right for every network are strongly connected and are the same since they are all 
generated by the reversible reactions 2F ^ G. By highlighting the nodes appearing in the same absorbing 
strong component, we see that in the state diagram of network 1, each fast component has a unique absorbing 
strong component. Comparing the top and center graphs in Figure we see the only difference occurs in the 
middle component. Network 2 has one more absorbing strong component than network 3, and the remainder are 
the same. 

We perform the reduction method on the three networks and obtain the reduced slow transition matrices for 
each network as follows. The detailed steps in the reduction are shown in supplemental material available at 
http://math.umn.edu/~othmer/Reduction.pdf, One finds that the columns of the 3 x 3 matrix 
Ki, the generator of the slow dynamics for network 1, are given by 


{ki)i 

- — feg fell/C5/C7(fe3+/C2 ) — fe8fe2fe5fe7 (felO+fell )—fe9 fell ^3^4 (fee+ ^ 7 ) 

(fel0+fell)(fe3fe4fe6+fe2fe5fe7 + fe3fe4fe7 + fe3fe5fe7) 

fegfeii 

feio+feii 

_ ksk2kgk7 _ 

fe3 fe4 fee fe2 fes fe7 kg fe4 fe7 “1“ kgkgkj 


{ki)2] 

_ fe8fe2fe5fe7(fe3fe5fe7 + 2fe2fe5fe7 + fe3fe4fe7 + fe3fe4fee) _ " 

fe2fe5^7(^3^4fee+^2fe5^7 + ^3^4fe7 + ^3^5^7)+|^i(/e4fee+^4fe7 + ^5^7)^ 

_— fe8fe2fe5fe7(fe3fe5fe7 + 2fe2fe5fe7 + fe3fe4fe7 + fe3fe4fee)_ 

fe2fe5fe7(fe3fe4fee+fe2fe5fe7 + fe3fe4fe7 + fe3fe5fe7)+^fei(fe4fee+fe4fe7 + fe5fe7)^ ’ 

0 


(ki)3 


■ 6fe9feii(feio+feii) “ 

3fefQ+6feioA:ii+fef^ 


0 


-6fe9feii(feio+feii) 

3fefQ+6feioA;ii+fef^ 
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12 ) 



Figure 5: The graphs of the full state transition diagram with initial state ®-(l, 0,0,0, 0, 2, 0) for network 1-3, red arrows denote 
slow reactions, (a) Network 1. (b) Network 2. (c) Network 3. Red and green arrows denote slow and fast reactions, respectively. 


The first three columns of K 2 are given by 


{K 2 )i 

-ih + h) - 


+ 


^3^2 

^2T^4 


(^ 2)2 

k^k2k7 


^3^4 

k 2 +kA 

^ 10+^11 

0 


/C 8 


{k2^kA){kQ^k7) 

—^5^2^? _ fegfeii 


{k2-\-k/C){kQ-\-k7) kiQ-\-kii 

0 

kio+kii 

0 


(^ 2)3 

2 k 8 

0 

■(^3 + 2 /^ 8 ) + 

^3^4 

^2 T^4 

0 


k^k 2 

k2~\~k4 
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Figure 6: The three fast components for network 1 (top), network 2 (center) and network 3 (bottom). The first 
has three fast components and three absorbing strong components, the second has three fast components and six 
absorbing strong components, and the third has three fast components and five absorbing strong components, 
and the last three columns are given by 

{K2)a 
0 
kg 

_ k^k2kr _ 

{k2+k4){kQ-\-k7) 

{k2 + k^){kQ+k7) -h ^2 + ^4 

k2+kA 

0 


—2k^k7 

k6+k7 


iK2h 
0 
0 
0 

2^5 ^ 2^7 

{k2+k4){k6+k7) 


+ 


2^5^4^7 

{k2+k4) {kQ-\-k7) 

0 


(i42)6 

6fc9fcii(fcio+fcii) 

sfe^Q+efeiofeii+Zcf]^ 

0 

0 

0 

0 

—6fc9fcii(/eio+fcii) 

3fcfQ+6fciofcii+/cfi 
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Finally, 



kgkii 

0 

2 ks 

0 

6 A; 9 A;ii(A:io+A^ii) 


kio^kii 

3 fc^Q+ 6 fcioA:ii+fcii 


0 

kgkii 

kio+kii 

0 

ks 

0 

K 3 = 

kgkii 

kio+kii 

0 

- 2 ks 

0 

0 


0 

kgkii 

kio+kii 

0 

-ks 

0 


ks 

0 

0 

0 

efcgfcii (fcio+fcll) 
3 fc^Q+ 6 fcioA:ii+fc^^ 

A number of general facts emerge from this example. 





1 . If each fast component has a unique absorbing strong component, entries of the left eigenvector of the fast 
transition matrix are all one’s , i.e. L = [ 1 , 1 ,...,!]. In this case, there is no need to calculate A's 
and S’s. 


2 . The column sum of Li is always one. 

3 . The dimension of the reduced slow transition matrix K is the total number of absorbing strong compo¬ 

nents. Our reduction method does not require the uniqueness of absorbing strong component in each fast 
component. Without loss of generality, we consider the middle component of network 3 for example, 
see Figure]^ It is obvious that the probability of ® choosing to transit to @ is fc2/(/c2 + k^), and the 
probability of transiting to @ is k^/{k2 + k^). Conceptually we can suppose that this 5 -state component 
is comprised of @ ^ @ ^ ® with probability k2/{k2 + k^), and occurs as @ ^ (Z) ^ ^ ^ ® with 

probability /c4/(/c2 + /C4), and the reduction method reflects this correctly. 

To see this, denote the probability vector corresponding to the 5 -state component by 

P=[Pn,Pl4,Pl9,P7,P3f- 


Then the left eigenvector L is found to be 


1 0 0 
0 1 1 


J2 


k 2 


^2 + ^4 k2-\-k4 

k^ k 4 


k2~\~k4 k2~\~k4 

and the probability vector of the reduced network is given by 


P^LP^ 


Pn + -i;^(Pr + Ps) 


k 2 +k 4 


Pu + Pn + A:{P7 + P3i 


k 2 ~\~k 4 


the first row of P corresponds to the probability of being trapped in the absorbing strong component 
whereas the second row corresponds to the probability of evolving to the other absorbing strong 
component The accumulated state associated with @ ^ @ ^ ® is obtained by switching off 

the reaction with rate ^4, which is (2, 0 , 0 ), and similarly @ ^ (Z) ^ ^ ^ ® gives rise to (2, 1 , 0 ). 


4 . Comparing the reaction network 2 and 3 as shown in Figure]^ we can see network 2 has two more reverse 
slow reactions, i.e. reactions with rate ks and k^. Consider the reduced slow transition matrices K2 and 
K3, if we let ks and k^ be zeros in K2, then K2 and are the same except K2 has one extra row and 
column of zeros, which correspond to the absorbing strong component that only exists 

in network 2. 
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Figure 7: Illustration of multiple absorbing strong components within one fast component: Network 2 middle graph 


4.5 The mean and variance 


Often the first and second moments of the distribution of each species is desired, and from the previous analysis 
one can find an approximate expression for the mean and variance of each species as follows. Define a projection 
operator Tk as 

Tk{n) = Uk, 

where n = (ni,..., Ug). Observing from the previous analysis that each state Sij, the state in Ci, has a 
probability distribution 7TijPi{t) at time one can find for each A: = 1,..., 5 

E[nk,t] = 'Y^Tk{'E^)pi{t)'Kij, (44) 

i,j 

Var[nk,t] = Y^[Tk{n'^)]‘^Pi{t)7rij - E[nk,t]^, (45) 

hj 


where denotes the vector n of molecular numbers corresponding to the state Sij. Thus we can obtain formal 
solutions ( |4^ and for mean and variance of each species on the slow dynamics. 

In linear models, it is possible to obtain explicit expressions for mean and variance [Gadgil et al, 2005| . In 
that case one can consider transitions between species as random walks of a molecule. This fact leads to a simple 
stochastic algorithm for computation and especially, if each fast component Ci in the reaction network is strongly 
connected, the quasi-steady-state distribution of the fast component is multinomial [Gadgil et al., 20051 . Thus 
one can compute the conditional expectation E[7l^{n)\h] as follows; The quasi-steady-state probability vector 
TTi of a component Ci can be computed by using 


KfTTi = 0,y^7r, 




= 1 , 


where ivij is a steady-state probability of species Mij in the component. If the reactant of a slow reaction 
is Mij, then one can obtain the transition rate for the reaction 


E[JZ^{n) I h] = K^TTijfi^ 

where kg is the transition rate constant of the reaction R^. 

Further insight into the structure of the slow dynamics is gotten from the fact that entry of is a 

conditional expectation when all fast components are strongly connected, as proven next. 
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Theorem 14 Suppose that fast components are strongly connected. Then the approximate transition rate TZi 
for the transition rate of all slow reactions from to fast component is 

Hi = 'y2iCiE[hi{n) I n] = ^£'[7^|(n) | n], 

where q is the transition rate constant of the slow reaction which transits states from to i^^fast component 

and hi is the combinatorial number of reactants of the slow reaction. 

Proof 4 First note that by an occurrence of the slow reaction, a state in the fast component £ ^ is transformed 
into a state in the fast component When a fast component is strong the variable h — A fU is uniquely defined 
in each component, we suppose that h — h} for Dj and h — h^ for Let be the vector of molecular 
numbers that denote the state Sik, the state in the component jCy 

Since LK^H is a Markov chain generator, it follows from the reduced equation ( |j9| ) that the approximate 
transition rate for a transitions from the to the fast component is given by the entry of LK^H, 

which is computed as follows. 

{LK^U)ij = {diag[Li ,..., diagfli,II2, • • •, ^nf])ij 


[ 11 ... 1 ] 

■ •• 

(^ 4)21 •• 

• {KtAm, ■ 
{KfY)2mj 


TTji 

7rj2 


AKtj)ma •• 



_ ^jrrij _ 


= 'Fj = Y1 cehi{n^’‘^)P{n = \ h = h^) 

q p q 

= ciE[hi{n) I h — fiY — E\TZl{ri) \ h — hP] 

where the block matrix of K^, Kfj, is a matrix of transition rates from the states of to those of 

and so [K^Yjpq is the rate of transitions from the state in into the state in jCy In the fourth equality, 
we made use of the fact that JYp{^ij)pq — JYi Oihi{n^^^) and iVjq is the conditional probability P{n = | 

h = h^). 

Remark 15 The result is true in general but the proof is somewhat more involved since the structure of the Li 
and the Hj is more complicated. 

When h defined earlier is uniquely defined in each discrete simplex we have the relationship between 
the full and reduced systems shown in the following table. 



Original system 

Reduced system 

variable 

n 

II 

stoichiometry 


iJs" = 

transition rate 

7^"(n) 

n^n) = 
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Moreover, under the assumptions in Theorem[^we can obtain an approximate chemical master equation on the 
slow time scale when the fast components are strongly connected. If the original master equation is given by 

jP{n,t) = ■ P{n-v£j^i,yt)-n{{n) ■ P{n,t)] 

+ IZ [^k{n - i^Elk)) ■ Pin - vElkyt) - niiri) ■ P(n,t)], 

k 

then the master equation in the reduced system is approximated (with error less than 0(e)) 

= y] P|(n - (fc ))pih - v£^ (fc) ,t)-Rlin)p{h, t) , (46) 

k 

where h — A^n, ^^(k) ~ ^ section we obtain a modified 

stochastic simulation algorithm from ( |^ . If there are sources or internal strong components the master equation 
for the slow system contains additional terms. 

5 A stochastic simulation algorithm for the slow dynamics 

5.1 A stochastic simulation algorithm based on the QSS approximation 

Preparation: 

• Identify fast and slow reactions. 

• For a given initial state n(0), switch off the slow reactions and generate the transition matrix of the 
fast component that n(0) lies in by implementing Algorithm[^ 

• Identify the fast simplex that n(0) lies in to determine 

• Define h = = A^^^^vE^. 

• Denote the fast simplex that n(0) lies in to be n(0) = A^^^^n{0) = hj, i.e., jth fast simplex and its target 
fast simplexes through slow reactions, denoted as (could be more than one). 

• Identify the fast transition matrices within the ith and jth fast simplexes, respectively, denoted as . 

Identify slow transition matrices from jth to ith fast simplexes, denoted as 

• Compute Li and Ilj such that LiK ^= 0 and = 0. 

Step 1. Scheme for a fast simplex : Computing slow reaction rates from the QSS state of fast a simplex 

1. Compute the slow transition rate iZ^{hj) = K^j = LiKf -flj from the jth fast simplex to all its target 
simplexes. 

2. Compute P|(nj). 

Step 2. Scheme for simulation of slow reactions : Simulating the slow reactions with the reaction rates 
obtained in Step 1 

1. Generate two random numbers ri and r 2 from the uniform distribution on (0,1). 
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2. SetT = 




log(ri) 


R: 


and choose k such that 


Tin 


tot 


Step 3. Update 

Let t ^ t + T and h ^ h + v£^. Go to Step 1. 


Remark 16 

7. The main idea of the modified Gillespie algorithm is to perform the traditional algorithm on a reduced rather 
than the original reaction network. For instance, in Example\W[ 


Af 


1 1 0 
0 0 1 


n(0) = A^n{0) = (2,0) o£^ = A^oS^ 


-1 1 
1 -1 


and u£^ provides the stoichiometric matrix for a reduced reaction network. The reduced reaction network can 
be written as 

Ml ^ M2’> 

where Mi, M 2 cire what we call “pseudo species” in the reduced network, and the new initial state n(0) = 
(2, 0) corresponds to only two molecules of species Mi. Each pseudo species represents one fast (absorbing) 
strong component, in this case 

Mi:A^B M2:C. 

The slow reaction rates in the reduced network are computed using LK^H as before except and are for 
the original reaction diagram rather than state diagram, i.e. 



A 

B 





c 

-h 

k2 ’ 

= C 

0 

h 

-k2 




Kii = C 


A B 

ke ks 


C 


^ 1,2 


A 

B 


k^ 

k^ 


After finding the forward reaction rate from Mi to M 2 to be {kik^ + A^ 2 ^ 6)/(^4 + ^ 5 ) the backward rate 
/c 4 + k^, one can initialize the Gillespie algorithm as usual. 


Example [T0| is the special case in which each fast component is strongly connected, i.e. only fast absorbing 
strong components or sinks exist. In general, one should identify sink, internal, and source strong components in 
the reaction graph. Once they are identified, the graph becomes a tree. Then one can use for example depth-first 
search to identify the fast discrete simplexes as defined inj^in order to find A^. 

The algorithm entails computation of the matrix of the slow dynamics, and whether one does this once for all 
initially, which may be advantageous when doing multiple realizations, or on the fly when only a few realizations 
are desired, is a matter of choice. In single realizations the entire state space generated by Algorithm 1 may not 
be explored and ’on-the-fly’ computation of the slow transition rates may be advantageous. 


6 Applications 

In this section we analyze three examples: a bacterial moter model, an enzyme-inhibitor model and a model 
from PFK system. These examples illustrate the reduction when there are relatively few states in the system. In 
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these cases one can find the approximate probability distribution of the slow variables by solving the reduced 
matrix equation dp/dt = LK^lip directly. Of course the main step is to find the matrix in the reduced equation 

I= 

Example 1 (A bacterial moter model) This example arises as a model for control of the rotational bias in the 
flagellar motor o/E. coli HOthmer, 2005]! . At the base of the motor are sites at which the protein CheYp can 
bind, and the occupancy of the sites biases the probability of switching the direction of rotation of the motor. 
Here we consider the following scheme for binding CheYp (represented byY)to these sites. 



nkiY 


CWo 

k-i 

CWi 

oq 4-t do 

nksY 

Oil tt di 

CCWo 

k-3 

CCWi 


{n-l)kiY 

2k-i 

(n 

2fc-3 


CW 2 
02 4-t P 2 
CCW 2 


2kiY 

(n — l)k-i 

2k^ 

(n - l)fc_3 


CWn-l 
O-n—l dn—1 

CCWn-1 


kiY 

CWn 

nk-i 

an tt/3) 


CCWn 

nk-3 



Here CWk and CCWk represent clockwise and counterclockwise flagellar rotation, respectively. Yle assume 
that the horizontal transitions are fast, while the vertical transitions are slow, which leads to two strongly 
connected fast components comprising the horizontal steps. 

Thus the two (n + 1) x (n + 1) fast reaction rate matrices are given by 


—nkiY k-i 0 


nkiY -{k-i + {n-l)kiY) 


{n-l)kiY 

II 

nk-i 


1 - 

1 

1 



Ki = 


Moreover, the slow reaction rate matrix is 


-nkzY k-3 

nkzY -(k -3 + (n - l)k 3 Y) 
(n - 1 )^ 3 ^ 


0 


kzY 


= 


-A B 
A -B 


where 


A = diag{ao,..., an) and B = diag{j3o,..., /3„). 

Let n = [Hi I 112 ] and let Hi, i = 1,2 be stationary distributions of xf, i = 1,2 respectively , i.e, 11, is the 
right eigenvector of xf corresponding to the zero eigenvalue with Iljj = 1, where Iljj denotes the entry 
of the vector 11,. Then the reduced equation can be written 

I--'"- 


where 


LX^Yl = 


-Er=i Er=i A-in2i 

E-=>i-inH -E-=i'A-in2i 


nk-3 

-nk -3 
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Thus we reduce the system into two-state model, 


ccw ^ cw, 

k~ 

where k^ = k~ = Here we note that k^ are functions ofki^ n, a, /3 and 

Y. These parameters are reported in the experimental literature. 

Figure^illustrate the evolution of probabilities of CW and CCW obtained from the full and the reduced 
equations. 






Figure 8: Shown here are the time evolutions of the probability that the motor rotates clockwise {Pr{CW) = 1) and counterclockwise 
{Pr{CCW) — 1) with different initial conditions using/ci = 2,A;_i = l,/c 3 = 2,/c _3 — ai — 0.1, = 0.08, F = 100. In 

the upper two figures, the initial condition is that the motor is rotating clockwise with no Y bound. In the lower two figures, the 
initial condition is that the motor is rotating counterclockwise with no Y bound. Blue lines are generated by computing the full system 
{Pr(CW) — Pr(CWi), Pr(CCW) — Yl Pr(CCWi)), whereas red circles are for the reduced system. 


Example 2 (An enzyme inhibition model) We consider an enzyme-substrate reaction system with a competi¬ 
tive inhibition 


ki /C4 

E + SZ^ES%E + P, E + I^EI, (47) 

k2 h 
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where I and P denote enzyme, substrate, inhibitor and product, respectively. Let ni{t)^i — 1,..., 6 denote 

the number of molecules of ES^ /, El and P, respectively. We assume that the two reversible reactions 
E -\- S ^ ^ ES.^ E -\-1 ^ ^ El are much faster than the irreversible reaction ES^E + P. After finding 


1 0 
0 0 
1 0 
0 1 

one can identify the slow variable 

n = A^n = (ni + ns + 77 . 5 , 02 + ns, 77.4 + ns, uq)^ . 

Since the fast subsystem E -\- S ^ ^ ES^ E -\-1 ^ ^ El has a deficiency of zero and is weakly reversible, one 
can use the result by IjAnderson & Kurtz, 2011^ for finding the equilibrium probability of the fast subsystem; For 
the convenience of computation, we assume that ki = ^2 and ks = ^ 4 . If the deterministic equilibrium values 
of E^ ES^ I and El for the fast dynamics are denoted by ci, C 2 , cs, C 4 and C 5 , respectively, one can find the 
only one solution 


10 10 
0 110 
0 0 0 1 
0 0 0 0 


and 


C 3 


C4 


0^2(1 + 0^1 + 0^2 + + (1 — 0^1 + 0^2 + ^1^3)^) 

2(0^2 + <^3) 

0^3(1 + Oi\ + 0^2 ~\~ 0^3 ~\~ \/ 4 q^i + (1 — 0^1 + OL2 + <^3)^) 

2 ( 0^2 + <^ 3 ) 


Cl = 0^1 — C3 — C5, C 2 = 0^2 — C3, C4 = CI3 — C5, 


where ai = ci + C 3 + C 5 , 0^2 = C 2 + C 3 , and 0^3 = C 4 + C 5 are conserved quantities in the fast subsystem. Using 
the result of HAnderson et al., 2010j , one finds the equilibrium probability of the fast subsystem as 


5 

p(ni,...,n5) = 

1=1 



n^ = 0 , 1 , 2 ,... 


where M is the normalizing constant for ~ ^ • • • 5 '^5 satisfy the conserved quantities ni + 

^3 + ^ 5 . '^2 + n 3 and + ns. The following figure shows the simulation results for the slow variable (the 
number of product molecules) obtained from the approximate algorithm and the exact algorithm. 


Example 3 (A model for the PFK system) We consider a reaction network from a model of the PFK step in 
glycolysis HOthmer & Aldridge, 1978]! 


Ai + Ei 
Ai + El 
A 2 + E 2 


k-i 

k-3 

k —s 


EiAi '%Ei + A2 
E*Ai 11 E^ + A 2 
E 2 A 2 E 2 + Product. 
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Figure 9: Enzyme-substrate model with an inhibitor. Comparison of approximate stochastic simulation al¬ 
gorithm (red dotted) to exact stochastic simulation algorithm (blue solid). Evolution of means and standard 
deviations of numbers of the product P when the initial condition (E’, 5, ES^ /, EI^ P) = (5,100,0, 5,0,0) 
and ki — k 2 — k 4 ^ — — 10, k^ = 0.1. The results are based on 5000 realizations. Concerning the relative 

CPU time for one realization of the stochastic simulation, if the approximate algorithm takes 1 second, the exact 
algorithm takes about 3.6 seconds, when a quad-core machine with Windows 8.1 and MATLAB 2014 is used. 


Here Ai^ A 2 , Ei^ E^ and E 2 denote EQP^ ADP, the low activity and activated forms of free PFK and the 
enzyme for the ADP sink reaction, respectively, and EiAi, E^Ai and E 2 A 2 represent enzyme-substrate com¬ 
plexes. 

If we assume three binding/unbinding reactions are much faster than others, there are two different fast 
subsystems Ci and C 2 in the reaction network. 

Cl : Ai + Ei^EiAi, Ai + El^ElAi 

k-i k-s 


and 


k^ 

C2 : A2 + E2 i — E2A2. 

k-5 

First we can find stationary distribution ofC 2 using the hypergeometric functions (fidcQuarrie, 1967]/ . Let initial 
numbers of A 2 and E 2 be uq and h^, respectively and define Q — Ifh{)> uq, then the stationary marginal 
distribution of A 2 is given by 


p (^0 + <^o)(Q'0 + cq — 1 ) • • • (uq + cq — A: -h 1 ) 

k\ (fio — ao + l)(^o — uo + 2) • • • (60 — clq + k) 

where A: = 1 ,... , ao + cq and D is a normalization constant. 

Similarly, if oq > bo, then the stationary marginal distribution of E 2 

p (j^) _ Co) (60 -h Co — 1 ) • • • (60 + Co — A: -h 1 ) 

^ k\ (ao — bo -\- l)(ao — 60 + 2) • • • (ao — bo A k)^ 


where A: = 1,..., 60 + co 

Since the fast subsystem Ci has a deficiency of zero and is weakly reversible, we can find the equilibrium 
probability similar to the enzyme inhibitor model; If we assume that ki = k 2 and ks = k/^ and the deterministic 
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equilibrium values of Ai^Ei^EiAi^El^ E^Ai are denoted by Ci A — 1,5, respectively, then we can find the 
equilibrium probability of the fast subsystem Ci 

p{ni,.. . ,n 5 ) = mJ]^ -^,ni = 0, 1 , 2 ,... 

*=i 

where riifi = 1,5 are the number of Ai^Ei^EiAi^ E^ E^Ai, subject to the conserved quantities ni+n^ + 
^5. '^2 + ^3 77.4 + 77.5, M is the normalizing constant for ' ff ^ p { n ) — 1 and 77.1,77.5 satisfy the conserved 

quantities 77.1 + 77.3 + 77.5, 77.2 + 77.3 and 124 + 77.5. 

0^2(1 + 0^1 + 0^2 + <^3 + ^ 4 o^i + (1 — 0^1 + 0^2 + <^3)^) 

2(0^2 + <^3) ’ 

0^3(1 + 0^1 + 0^2 + <^3 + ^ 4 q^i + (1 — 0^1 + 0^2 + <^3)^) 

2(0^2+ <^3) 

Cl = — C3 — C5, C 2 = 0^2 — C3, C4 = CI3 — C5, 

= Cl + C3 + C5, a2 = C2 + C3 and 0^3 = C4 + C5 arc conserved quantities in the fast subsystem. 

The following figure shows the simulation results for the number of product, which is a slow variable. 




Figure 10: A comparison of the approximate stochastic simulation algorithm (red dotted) to the exact stochastic 
simulation algorithm (blue solid) for the model of PFK reaction system. Evolution of means and standard 
deviations of numbers of the product when the initial condition (Ai, E’l, E’lAi, E’^Ai, A 2 , £' 2 , £ 2 ^ 2 , P) — 
(100, 5,0, 5,0,100, 5,0,0) and k 2 = = kQ = 0.1, ki = k-i = ks = k-^ = k^ = k-^ = 10. The results are 

based on 5000 realizations. Concerning the relative CPU time for one realization of the stochastic simulation, 
if the approximate algorithm takes 1 second, the exact algorithm takes about 2.9 seconds, when a quad-core 
machine with Windows 8.1 and MATLAB 2014 is used. 


7 Conclusion 

We developed a reduction method for stochastic biochemical reaction networks with coupled fast and slow 
reactions, and formulated an associated extension of the Gillespie method. The reduced equation for the time- 
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dependent probability distribution on the slow time scale involves a projection in state space, and we show that 
the generator of the reduced system is also Markovian. We found that the transition rate for a reaction on the 
slow time scale is approximated by the expectation of the original transition rate conditional on the invariant dis¬ 
tribution of the fast dynamics. Throughout the reduction the invariants of fast subsystems play a significant role, 
similar to what was found earlier in the analogous reduction of deterministic systems [Lee & Othmer, 2009 1. 
When the stationary distribution of the fast dynamics can be computed, these can be used directly in a modified 
stochastic simulation algorithm. Several biological examples, including a model of motor behavior for a single 
flagellum, an enzyme-inhibitor model and a model for the PFK system, illustrate the numerical accuracy of the 
approximation for two lowest moments. 

8 Appendix 

8.1 The general case in Example [l0| 

When there are Nq molecules, there are A^o +1 fast components and the transition matrix K is block tridiagonal, 
i.e. 


K = 


1 


k( + Kf K 


^ 2,1 


1,2 


-Ki + KS 


^3,2 


where the fast blocks are given by 


^2,3 


K 


Nq+I.Nq 


-kI 

e 






K( = 


-Noki fe 

Noki -[{No - l)ki + k 2 ] 2k2 

{No - l)ki -[{No - 2)ki -h 2fe] 3fe 


Nok2 


ki -Nok2 


(Aro + l)x(iVo + l) 


N _ 


-(iVo - l)ki k 2 

{No-l)ki -[{No - 2)ki + k 2 ] 2fe 

{No-2)ki -[{No - 3)ki + 2k2] 3k2 




—ki k2 
ki —k2 


.^lo+l=0- 


(iVo-l)fc 2 
fcl -{No - l)fc 2 


NqXNq 
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The slow off-diagonal blocks are given by 


■ fcs 

/C4 k^ 

II 

"3 ri 

’ 2/C5 

2/C4 2/C5 

, • • • ,^^o,A^O + l — 

iVofe 

1 - 

1 _ 

(A/'o + l) X A^o 

1 

to to 

1 _ 

Arox(iVo-i) 

NokA 


For the lower diagonal blocks, 



Noke ks 


(Nq — l)ke ks 

Kli = 

{No - l)k6 2k3 

KI 2 = 

(No - 2) fee 2k3 


ke Noks 

Nq X (A?'o + l) 

ke {No - l)k3 


^^ 0 + 1 ,[ ^6 ks Finally, the K? along the diagonal are diagonal matrices of the same dimension 
as Ff/, and the (j, j)-th entry of is the negative sum of rates leaving j-th node of the i-th fast component. 

In this case the transition rate is given by 

Kj - 


8.2 Moment equations of the invariant distributions 

The low-order moments of the distributions for the fast systems play a role in the QSS reduction in Section]^ 
and here we consider the low-order moment equations. 

Theorem 17 Let r be the total number of the reactions in the system. Then the invariant (steady-state) distribu¬ 
tion of P{n^ f), which we denote P{n), satisfies 

r r 

ni{n - iyS(^i)P{n - TZe{n)P{n) ( 48 ) 

e=i e=i 

and and the first two moment equations lead to 

EfiSnin)] = 0 , ( 49 ) 

r 

\y£{i) ® E[riJZi{n)] + E[riR.i{n)\ ® ® oSf^i^ElTlfin)]] = 0 . 

1=1 

(50) 

Proof 5 At the steady-state 

r r 

y^7e£(n- i/f(^))P(n- i/f(£)) = ^ni>{n)P{n) (51) 

t=i t=i 

By multiplying by n and summing over all the values ofn G we obtain 

r r 

EE nUfin - o£(^())P{n - o£(^i)) = EE nPi{n)P{n). 

n i=l n 1=1 
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Using the transformation n — 08 (^ 1 ^ n on the left side, we obtain 

r r 

+ i'S^i))ne{n)P{n) = nTZein)Pin). 

n 1=1 n i=l 

By subtracting the right side from the left one, 

r r 

1=1 n i=l 

Thus we conclude that 

v£E\R{n)\ = 0. 

If the deficiency S = p(£) — p{o£) is zero, then E[JZ{n)\ is a cycle in the graph fjOthmer, 7979|/ j^ 

At the next order we multiply equation by a tensor product n ® n and sum over n. Then by a similar 
argument, we obtain 

r 

® ^ ^ 77 / 7^2 (^)^(^) ^ ^ p(^ifj ® “ 1 “ ® ^ ^ p 

i=l n n n 

r 

^ E ® ElnUiin)] + E[nTZi{n)] ® (g) u£(^i)E[TZi{n)]] 

1=1 

= 0. 


r 

E 

i=l 


If 5 = 0, the two lowest moment equations can be simplified to 

E[n{n)] = 0 (52) 

\y£[i) ® EfiEfin)] + E[nEi{ri)\ ® = 0. (53) 

When all reactions are linear the problem is much simpler, and the evolution equations for the first and second 
moments can be written explicitly in terms of those moments [ Gadgil et aL, 2005 1. 

As a consequence of Theorem [TT} similar equations can be obtained for the quasi-steady-state of the proba¬ 
bility distribution for the fast subsystem in a two-time scale stochastic network. We first define the expectation 
of a function f{n) over a discrete reaction simplex Cf for the fast subsystem as follows. 

ECf[f{n)] = 

neCf 


Corollary 18 Let Vf be the total number of fast reactions. Then at the steady-state of the fast subsystem, the 
governing equation is given by 

y] n{ (n - u£(.^)Pin - v£(^^) = n{{n)P{n) (54) 

1=1 1=1 

and for each discrete reaction simplex Cf, 

Ecf[npn)] = Q, (55) 

1=1 

^The reader can show that the presence of inputs or outputs of the form given in Tablej^does not alter the deficiency. 
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Proof 6 The 'state-wise' form of the master equation ([79|) can be written 


^P(n, i) = XI ^ ~ '^e (”) ’ ^)] 

+ X [mn - • P{n - vS^yt) - nUn) • P(n, t)], (57) 

k 

where TZ^ and TZ^ are the transition rates of fast and slow reactions, respectively and and £^ are incidence 
matrices for fast and slow reactions, respectively. 

In the previous theorem, substitute £^ ^ 7Z^ and [•] into £^ TZ and E[-] and use the full rank assumption 
on . 
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